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^ ; Abstract 

The two-dimensional extended Hubbard model that includes a nearest- 
neighbor Heisenberg interaction is studied using a mean-field theory where 
quasiparticles are defined by an U(8) group of canonical transformations. 
The theory is a generalization of the ordinary BCS theory, and Balian and 
Werthamer's theory of 3 He that permits both broken gauge, spin and sub- 
lattice symmetry. This allows us to investigate superconductivity, antiferro- 
magnetic order, charge density waves and, by twisting the spin quantization 
axis, spiral antiferromagnetic order in the same theory. Our results for pos- 
itive Hubbard U and Heisenberg exchange J suggest that antiferromagnetic 
ordering dominates close to half filling, while spiral states and d-wave super- 
conductivity compete when doping is introduced. For moderate values of J, 
we find a phase diagram where a phase transition occurs from an antiferro- 
magnet to a d-wave superconductor as doping is increased. A narrow region 
of (s + id)-wave superconductor is found for some values of J and U. 

PACS numbers: 74.20.Fg, 74.25.Dw 
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I. INTRODUCTION 

For many years the two-dimensional (2D) Hubbard and "extended" Hubbard model have 
served as simple models of high-T c materials and a paradigm of strongly correlated electrons. 
However, despite many creative attempts to develop approximations to attack the problem, 
the 2D Hubbard model has defied a definitive analysis. 

In contrast to many of those efforts, the goal in the present paper is not to try to model 
the behavior of high-T c superconductors or even to make particularly strong statements 
about the 2D Hubbard model. Rather, we wish to make a definitive analysis of the BCS 
mean-field theory, an approximation method that has been remarkably successful in gaining 
information about simpler model systems with both broken spin and gauge symmetry. It 
is clearly important that we know what this simple approximation has to say about the 
Hubbard model before we can be confident in applying more sophisticated techniques. 

Although the program sounds straightforward, the method becomes surprisingly compli- 
cated if one insists on preserving the symmetries known to exist in the Hubbard model at 
half filling. The complication is due to the exact SO (4) « SU(2) x SU(2) symmetrytl pre- 
cisely at half filling and that this symmetry mixes many phases that are typically ignored 
in simpler calculations. To do a proper job, we must therefore include all the order parame- 
ters that have been discussed before for this system, among them spiral spin-waves,aH Neel 
antiferromagnetism,Q'S <i-wave and s-wave super conduct ivityoBi!3 as well as all other phases 
that are related to these via the symmetry group B 

A consequence is that a multitude of possible order parameters must be retained to 
provide a self-consistent theory. Our mean-field theory systematically enumerates a very 
large set of representations of the possible broken symmetries and, among other phases, 
allows for the possibility of all types of ordered phases that have been discussed in previous 
analyses, including antiferromagnetism, spiral spin-waves and <i-wave superconductivity. 

The mean-field theory that we employ is based on an U(8) set of canonical transforma- 
tions that are analogous to the U(2) canonical transformations discovered by Bogoliubov and 
Valatin in their analysis of the broken gauge symmetry in BCS theory.EJo Their approach 
was generalized by Balian and Werthamer to a set of U(4) transformations in their study 
of superfluid 3 He, where also spin rotational symmetry is broken.!^ We, in turn, extend 
their analysis to include broken sublattice symmetry, which extends the relevant group of 
canonical transformations to U(8). 

The organization of the rest of this paper is the following: After providing a brief review 
of related results in Sec. [TI], we analyze the symmetry of the Hubbard model in Sec. |TTT| . The 
representations of the SU(2) x SU(2) symmetry group is made via an 8 x 8 Clifford algebra, 
which provides a generalization of the Dirac and Pauli matrices used previously in the Nambu 
formalism applied in the study of BCS theory and 3 He.t3 In Sec. |V|, we derive the mean-field 
theory for the Hubbard model at half filling and introduce the point-group symmetry, while 
we in Sec. [V] extend the analysis to include the possibility of spiral antiferromagnetic states 
that have been proposed as likely phases away from half filling. We discuss the reformulation 
of the self-consistent equations in Sec. [VT|, and the numerical solution procedure in Sec. |V11 . 



Our results are presented in form of phase diagrams in Sec. |VIII| , before the final discussion 



in Sec. [TX| Appendices provide additional mathematical details. 



II. BACKGROUND 

The Hubbard model is given by the Hamiltonian 

H = H - fxN + # Hu bb , (1) 



where 



H = ~t £ (4, CT c R > + h.c), (2) 

<R,R'),o- 

iV = 5> Ri(T , (3) 

R,(X 

^Hubb = ^£(™R,T - |)Ku - \) > ( 4 ) 



R 



and where 7i Rj0 - = c R(T c Ro .. The sums in the hopping terms are summed over all pairs of 
nearest neighbors on the 2D square lattice. The form of the Hubbard term is chosen so as 
to make the Hamiltonian particle-hole symmetric, and the system is half filled (one electron 
per site) when the chemical potential \x is set to zero. 

Several extensions to the Hubbard model have been studied. To have the possibility of 
mimicking the behavior of higher order terms which would be present in a more sophisticated 
weak-coupling expansion, and to physically incorporate the effect of spin-spin interactions, 
we extend our Hamiltonian by adding the nearest-neighbor Heisenberg interaction i^Heis 
given by 

-f^Heis = J 2-^ ^ R ' ^ R ' ' (") 

<R,R) 

where Sr = |( c R + c R . )er( c R ♦ c R . ) T and a is the vector of Pauli matrices. The Heisen- 
berg term breaks neither spin nor pseudospin symmetry, and thus fulfills the most important 
constraints on effective higher-order terms in a theory for the pure Hubbard model at half 
filling. 

In principle, we should also investigate the nearest- neighbor Coulomb interaction. Al- 
though it may be physically important, it violates pseudospin symmetry which is already 
broken by the chemical potential term. Since it adds nothing theoretically fundamental to 
the model, we chose to not include it in our analysis. 

An appropriate starting point is the pure Hubbard model at half filling. Here, the model 
is symmetric under the SU{2) x SU{2) group of global spin and "pseudospin" rotationsBEJcHl 
and the phase diagram is rather well understood. 

For positive values of U, the ground state is expected to be a 2D antiferromagnetic (AF) 
(7r,7r), or Neel, state. For large U/t, second-order perturbation theory in t/U yields an 
effective theory equivalent to the spin-| Heisenberg antiferromagnet. The order parameter 
is given by the staggered magnetization, which is a vector which breaks the global SU(2) 
spin and the discrete translational symmetry of the underlying model. 

For negative U, the groundstate is either an s-wave superconductor (SC) or a (n, n) 
charge density wave (CDW). The real and complex part of the SC order parameter together 
with the CDW order parameter form a triplet which transforms as a vector under SU(2) 



pseudospin rotations. The groundstate we call "mixed" SC-CDW since each of these phases 
are degenerate and related to each other by a continuous symmetry. 

The cases of positive and negative U are directly related via the "Shiba" transformation 
ZB3 This discrete transformation has the effect of changing the sign of U in the Hubbard 
model, and also mapping rotations in spin space into rotations in pseudospin space and 
vice versa. Furthermore, a state with broken pseudospin symmetry for U < (SC or CDW 
state), has a direct mapping to a state with broken spin symmetry (an AF state), and since 
SC already appears for an infinitesimal negative value of U, AF order must appear for an 
infinitesimal positive value. Hence, there is no Mott transition as a function of U/t, and the 
Neel order persists for all positive values of U. 

Away from half filling, we are less certain of what phases may occur. The pseudospin 
symmetry is broken, with only the U(l) gauge subgroup remaining. Spiral spin-waves, 
where the antiferromagnetic order parameter twists with a pitch along a symmetry axis, 
has been suggested as a likely candidate. This phase has been seen in several theoretical 
calculations and have been suggested to explain the incommensurate spin correlations seen 
experimentallyaia We study this class of states by imposing a twisted spin quantization axis 
in the generalized Hubbard Hamiltonian before it is analyzed by the Hartree-Fock-Bogoliubov 
met hod. a 



III. GROUP THEORETICAL FRAMEWORK 

A. Using multispinors to represent symmetries of the Hubbard model 

Let us first consider the noninteracting theory H = H q at jjl = defined for a square 
system with N sites, periodic boundary conditions, and lattice constant 1. This is simply 
diagonalized as H = J2k e k c k o- c k <n where the k-sum runs over N points in the first Brillouin 
zone (BZ). The single particle dispersion relation is e^ = —2t(cosk x + cosk y ). To make 
manifest the particle-hole symmetry, we define the vector Q = (jc, it) and the operators 

a L = C L , when e k < , 
b l,a = Ck+Q,<r > when e k+Q > . 

In terms of these operators, Hq is written as 

H = ]r"ek(4,a a k, CT - & k ,A,a 
k,er 

- a -k,a a -k,- + b -k,a b -k,a) , ( 7 ) 

where the summation denoted by J2" runs over the four times reduced Brillouin zone cor- 
responding to ek < and k y > (see Fig. [I]). Note the special order of the operators in 
Eq- ©• 

We next introduce the "Shiba" transformation Z which acts on the position space creation 
and destruction operators c\ a through the canonical transformation cj , i— ► (— l) r c r ,, c\^ \— > 

c r p where the factor (— l) r = e tQi ' r induces a change of sign on one sublattice. The Shiba 
transformation is hence a particle-hole transformation, together with a local change of gauge, 
which only acts on the spin-down operators. 



In reciprocal space, spin rotations and Z act naturally on the 8-component multispinor 

^k — ( a k,t> a k,|' ^k t' ^kh °-k f a -kh -KV -k,j.) which carries definite momentum k mod Q. 
In this basis, Z is represented by an idempotent matrix whose entries are all zero except 
Z ltl = Z 3)3 = Z 5)5 = Z 7 j = Z 2 $ = ^4,6 = Z 6A = Z 8j2 = 1, and which acts on \I/ k by 
matrix-vector multiplication, \& k i— > ^ k Z. It is well known that Z is an exact symmetry of 
H but changes the sign of the Hubbard term [/.E§E3iH3 It is also simple to verify that the 
nearest neighbor Heisenberg term Hn eis is invariant under Z. 

A pseudospin transformation R' is defined as R' = ZRZ where R is an ordinary global 
SU{2) spin rotation. Since [H,Z] = and [H,R]=0, we have [H,R'] = 0, and since Z does 
not commute with R, we see that the entire symmetry group of H is S = SU(2) x SU(2). It 
follows further that because \l/ k defines a representation of both R and Z, it also determines 
a representation of S. 



B. A basis spanning the space of Hermitian 8x8 matrices 

To represent the action of spin and pseudospin transformations on the basis \l/ k , we 
introduce a set of seven Hermitian 8x8 matrices Pa constructed in blocks from the ordinary 
Dirac matrices 7^ as shown in Table |. By explicit computation, it can be verified that Pa 
form a Clifford algebrac3, i.e., (3a(3b + PbPa = ^Qab-, where qab = except for g 00 = 1 and 
9nn = -1, n=l,...,6. 

The Pa matrices have been constructed so that (/5i,/3 2 ,/3 3 ) transform as a vector un- 
der spin rotations and ((3^, (3 5 , (3 6 ) transform as a vector under pseudospin rotations, while 
{(3&,f3$,f3§) and (/?i,/?2, ^3) transform as a scalar under spin and pseudospin rotations, re- 
spectively. Taking multiple products of the matrices Pa, one can construct a complete basis 
for the vector space of 8 x 8 Hermitian matrices (with real- valued coefficients). The Clifford 
algebra contains four independent spin/pseudospin scalars, and it follows that there are also 
four independent sets of spin and pseudospin vectors and spin®pseudospin tensors, making 
a total of 4 x (1 + 3 + 3 + 9) = 64 basis elements which together span the space of 8 x 8 
Hermitian matrices." 

Our goal is to construct a basis that has simple transformation properties under spin 
and pseudospin rotations, and with the above classification scheme in mind, we label the 
64 base matrices by B", where < fi, u, k < 3. The upper index k enumerates each of the 
four independent sets of matrices associated with each of the four scalars. Transformation 
properties under spin rotations are indexed by /x, where scalars carry the index \x = 0, while 
fj, = 1,2,3 represents the components of a spin vector. Similarly v identifies pseudospin 
scalars and vectors. For example, B\ x transforms as a scalar under spin rotations and as the 
first component of a pseudospin vector. The matrices B* are constructed so that they are 
not only Hermitian, but also have natural transformation properties under parity, sublattice 
exchange and time reversal, and so that the indices \x and v are simply interchanged under 
Z. 

To explicitly construct the basis matrices we introduce the four scalars. Three of them 

f ry \ 

are the identity matrix 11, p , and T = 1P0P1P2P2, = i I, where 75 indicates the 

pseudoscalar Dirac gamma matrix 75 = ^70717273; the fourth scalar matrix is defined as the 
product PqT. The four scalars are denoted T K with the following identifications: T = —11, 



Ti = /3o, T 2 = — r and T 3 = — iT/3 . To make subsequent formulas simple we also define 
Q = Q = 1? ty — ftj, and flj = /5 J+3 , where j = 1,2,3. Finally, introducing r = 1 and 
Tj = y/—T, the basis matrices are defined by 

B£„ = T K ^Cl u , for k = or 1 , 

B;, v = T^T^Clv , for K = 2 or 3 , (8) 

where < //, v < 3. The phase factors r M are chosen so as to make B* Hermitian. 

Letting the subscript "m" denote the collection of indices k, /i, u, we find that the matrices 
B m are orthonormal in the sense that 

Tr (B m B m ,) = &8 m>m i , (9a) 

B m B m = n . (9b) 

It follows that an arbitrary Hermitian 8x8 matrix M can be expanded as 

M = Y J oTB m , (10) 

m 

where 

a m = §Tr(£? m M) . (11) 



IV. WRITING THE MEAN-FIELD HAMILTONIAN IN TERMS OF 8 x 8 

MATRICES 

The interacting Hamiltonian contains not only one-particle but also two-particle terms. 
We use Wick's factorization^ to express the expectation values of the two-particle terms as 
sums of all products of one-particle terms: 

(O x 2 OzO A ) = (O x 2 )(OzOt) + (0,0,) (0 2 3 ) 

-(0 1 3 )(0 2 4 ), (12) 

where O; denotes any creation or destruction operator. 

After transforming H + i^Heis m Eqs. ([!]) and ([5]) to reciprocal space and then applying 
the factorization we find 

(H) = ^2(e k - M)(ck, CT c kiCT ) + — Y, <K k i - k 2 + k 3 - k 4 ) [US^S^ - -(27 kl _ k4 + 7^-^ 

k,cr k^ ,k2,k3,k4 



x 



( C ki, CT C k 2 ,cr)( C k 3 ,cr' C k4,cr') + ( C ki,o- C k 4 ,cr' ) ( C k 2 ,cr C k3,cr' ) ~~ ( C kx ,a C k 3 ,(j' ) ( C k 2 ,o- C k 4i cr' ) J 5 (13) 



where N is the number of lattice sites, 7 k = cos(k x ) + cos(k y ), and 8 indicate a ^-function 
modulo the reciprocal lattice. 

Let us now consider the one-particle expectation values. In ordinary mean-field analysis, 
all one-particle expectation values that carry non-zero momentum are assumed to be zero. 

6 



Similarly, in our formalism the natural set of expectation values are related to the elements 
of the 8x8 matrix of operators 



yi ® * k ) 



a k,T a k,T a k,T a k,| • • • 

(14) 



a k,| a k,T 



t 
k 

V J 



whose expectation values are all of the form (c^^c^^), (4 i>CTi cL k2i(J2 ), and (c ki>CTi c_ k2(J2 ), 
where ki = k2 or ki — k2 = ±Q with Q = (jr, -k). Thus we allow all possible expectation 
values of operators that carry momentum (0, 0) or (n, it), and those that carry spin and/or 
charge to be non-zero. Expectation values carrying momentum (n, n) correspond to stag- 
gered order, net charge represent superconductors, and net spin is associated with broken 
spin symmetry. 

To make clear the relation between the 64 individual operators in (\l/ k ® ^ k ) and the 
irreducible representations (irreps) of SU(2) x SU(2) we define a set of 64 operators by the 
following linear combinations 

<=~Tr[B m (^®* k )]. (15) 

We then use these operators to rewrite our mean-field Hamiltonian. 

We begin by rewriting the quadratic (one-particle) part of the Hamiltonian in matrix 
form. From Eq. (0) it is easily seen that H = Ek'^kTr^^t. g> \p k )) = 8£ k "e k (c4 ) k , 
where the basis matrix Bq is diagonal and has diagonal elements (1,1,— 1, — 1, — 1,-1,1, 1). 
The expectation value of the entire quadratic part of the Hamiltonian then becomes 

(H - »N) = S^'^iiaM-^ii^osh^- ( 16 ) 

hop fill 

The interaction (two-particle) part of (H) may be evaluated in a similar manner. To sim- 
plify the corresponding expression, we introduce some shorthand notation. Let (a" ) kk / = 
(( a /L)k)X(a^i/)k')- With this notation, the expectation value of the Hubbard interaction 
becomes^ 

\#Hubb) 



LOU y-^i/ . 3 . 2 
nr Z^ \ a 0i/kk' 

k ' k ' SC-CDW 


" ( a J0/kk' + ( a 0J/kk' 


AF SC'-fiU 



-£&>k£, (17) 

FM 

where the index % is summed from 1 to 3. The terms representing the mixed superconduct- 
ing and charge density wave (SC-CDW), antiferromagnetic (AF) and ferromagnetic (FM) 
ordering are underbraced. The third ("z") component of the SC'-fiU term is recognized as 
the filling while the x and y components represent the real and imaginary part of a kind 
of staggered superconductor (SC). This term together with the first are the analogues of 
the FM and AF states under the change of sign of U. We find numerically that the order 
parameter for SC is always zero for the interaction parameters that we consider. 




To make a more systematic investigation of all order parameters, we note that the Hub- 
bard model on the 2D square lattice has the point group symmetry C± v , which has four ID 
irreps {A\, A 2 , B% and B 2 ) and one 2D irrep (E) (see App. |A] ). The total symmetry is 
hence C iv x SU(2_) x SU(2), and it is possible to perform a complete classification of the 
order parameters^ 

Since the Hubbard interaction is on-site, the interaction coefficients are independent 
of k and all the prefactors of the Hubbard term belong to the A\ irrep. Another basis 
function for this representation is 7k = cos(fc x ) + cos(k y ), which occurs in the hopping and 
Heisenberg terms and in "extended s-wave" superconductivity. The Heisenberg term has also 
components in the B\ and the E representations. The B\ terms (odd under k x <->• k y ) have 
the k-dependence r^ = cos(k x ) — cos(k y ) and the two basis-functions of the E representation 
(odd under parity, k — > — k) are (k x = sin(k x ) and (k y = sm(k y ). The expectation value of 
the Heisenberg term reads, in its whole glory, 

+ i(7k7k'+%%')(-3(«oo)kk'-3 (aoi)kk' + (a^o)kk' +(<4j)lk') 

hop SC-CDW spin nem. 

+U(k x (K + (k y (k> )(-3{^ )kk' - 3(a° i )^ k , + (a° ) kk , + K>kkO- (18) 



From this form we see that the Heisenberg term renormalizes the hopping and also affects the 
antiferromagnetic and ferromagnetic ordering. Apart from these terms, the most interesting 
term is the o^-term representing superconducting and CDW ordering. This term is split 
into two irreps, A\ and B\ corresponding to 7k and 7?k- The A\ part gives rise to "extended 
s-wave superconductivity" in certain regions of parameter space. The B\ part introduces 
ci-wave superconductivity, and the ^-component of this part is recognized as an orbital 
antiferromagnet. We also recognize the B\ representation of the af term, which represents 
a spin nematic state with the order parameter (c kQ c k/3 ) = ir/kCa/rd, where d is a real vector. 
Finally, there is a term representing p-wave superconductivity and since it is odd under parity 
it belongs to the two-dimensional E representation. Some of the other terms may have been 
discussed in the literature but, since we have found them not to be energetically favored, we 
do not discuss them further. 

In addition, the full Hamiltonian gives a term that shifts the chemical potential, and 
the Hubbard and chemical-potential terms produce constants that shift the energy. These 
terms have no significance in our analysis and they are neglected in the following. 

In collecting all these terms, it is useful to write the expectation value of the Hamiltonian 
in the following generic form: 

(fl)=E\Ak)+ E 'UA(OW), (is) 

(,m,k l,m, k.k' 

where a; jm and b^ m are coefficients of the linear and quadratic terms, respectively, and 
/3 k = l,7k,^k, • • • are trigonometric prefactors from various irreps of the point group C± v . 
The index "1" labels the irreps of C^ v , while the index "m" labels the irreps of SU(2) x SU(2). 



V. INCORPORATING THE SPIRAL PHASES IN THE HAMILTONIAN 

The most important unanswered question is what happens when the electronic density 
deviates from half filling. A widely discussed scenario is that the antiferromagnetic state 
adjusts to incorporate the excess electrons and changes into a spiral phase — an antiferro- 
magnetic phase whose order parameter shows a spiral spatial distribution. Given that the 
spiral antiferromagnet is the energetically preferred phase in some other mean-field calcula- 
tions, we must clearly incorporate this scenario in our analysis. Unfortunately, spiral states 
cannot be directly generated by the Bogoliubov transformations we have used this far. We 
therefore generalize our method by imposing a spiral twist of the quantization axis on the 
Hamiltonian before the mean-field approximation is applied. 

To describe any spiral we need a 2D vector q in the plane, which controls the direction 
and pitch of the spiral spin wave, and a 3D unit vector f2 which defines the axis around 
which the spin is twisted. The spin-twisting canonical transformation is then generated by 
the operator 

(4 >T 4jexp[*(q.r)(fi. CT )](^ , (20) 

where a is the vector of Pauli matrices. This canonical transformation is easily written as 
a unitary transformation of the creation/annihilation operators, 



[cos(q-r)a + zsin(q-r)fJ-<r] f). (21) 





Since the Hamiltonian is invariant under global spin rotations the direction of f2 can be 
arbitrarily set to be along z. 

When applying this spiral spin transformation to H, the Hubbard interaction is invariant 
since it is spin-rotation invariant and on-site, but the hopping term transforms into Hq = 
Ek £ k(%+q,t + n k-ci,i), which, using our operator definitions, is written as 

(H«) = -2t ]T "8[cos(g n )cos(£; n )((c4)k> 

k,n=x,y 

+ sin(g n )sin(A; n )((c4) k )]. (22) 

The Heisenberg term transforms into H^ cis = JJ2(r,r') Sr ' Sr/ cos[2q • (R — R')] + 
2>Sr>Sr' sin 2 [q • (R — R')], where S^ denotes the ^-component of the local spin operator 
at site R. Using the same notation as in Eqs. (|TTD and (18) with the extension that i is an 



index that is summed over 1 and 2 only, the Heisenberg expectation value is 

(«3o)kk' - («3o)kk' + c Os(<?n)[(«!o)kk' - (cX^L'] 



(-^Hcis) — ~irr 2.^1 



N 

k,k' 
n=x,y 



+-cos(A;„)cos(A4) j - [1 + 2 cos(2q n )][(a 1 00 )l k , + (a^ k ,] + (c4) kk , + K\) kk , 
+[2cos(2g n )-l][(^ } kk ,+ (a^) kk ,] 



+ sin(/c n ) sin(A4) - [1 + 2cos(2g n )][Ko)kk' + {<i)L>} + (<o)L,> + (<4)L 



+ [2cos(2q n }-l}[{al )i k/ + (at i )i k ,} 



(23) 



The large number of terms in this expression is due to the broken spin-rotational and point- 
group symmetries. 

VI. SOLVING FOR THE STATE OF LOWEST FREE ENERGY 
A. Self-consistent equations at finite temperature 

In the spirit of standard BCS theory we introduce the reduced Hamiltonian 

H = £ "ai, m PlaZ + £ \m/£&aIM , (24) 

£,m,k l,m, k,k' 

which has the generic expectation value given in Eq. ([191) . We also define the mean- field 
order parameters (gap functions) 

Ai,™ = EX« ) • (25) 

k 

Using the assumption that the fluctuations in the operators a™ from their mean-field values 
are small, we substitute 

<=«-«}) + «> (26) 

into H in Eq. ([24]), drop terms quadratic in (a^— (c^k)) an d fi n d the mean-field Hamiltonian 

H mi = J2 "( a i,m + 2b litn A ltm )fi t a£ - Y, A lm- (27) 

Z,m,k l,m 

Aside from a constant that is unimportant in this discussion, H m { can be recast in the form 

if mf = ^V(k)K) a W 

k,a/3 

= §£V[/ i (k)(*j c ®* k )], (28) 

k 

where 

h(k) = J2( a i,m + 2b l , m &i,m)PLB m . (29) 

Ijn 

Introducing the matrix of expectation values 

/ k =<(*i®*k)>. (30) 

10 



it follows from Eqs. (|15|) and (^5]) that 

\ m = |E>kTr (B m f k ) . (31) 

k 

To evaluate /k we note that the mean-field Hamiltonian is bilinear in \l/ k and can be diago- 
nalized, 



^ mf = E^(k)(xI) Q (x k r, (32) 

a,k 

by the canonical transformation 

xl = U k *l , (33) 

where U k is an U(8) matrix. Standard arguments from statistical mechanics give in the 
diagonal case 

((Xk®X k k 7 ) = (l + ^ Q(k) )-X 7 , (34) 

(f3 is the inverse temperature, f3 = 1/fceT) which, when transformed back to the operators 
\l/k, results in 

/k=((^ k ®^ k )) = (l + e /3/l «)- 1 . (35) 

Equations (|29|), (]3T]), and (|35| ) constitute the self-consistent equations to be solved for A/ :m . 
The solutions extremize the free energy F, 

F=(H}-TS } (36) 

where (H) now includes the chemical potential, and where 5* is the entropy. If more than 
one solution is found, then the one with the minimal value of F is the physical one. To 
calculate F, both terms in Eq. (|36"D must be explicitly evaluated. The evaluation of (H) is 
straightforward — using Eq. (|T9|) we see that it is equal to 



(H) = EK»A m + b lym AU . (37) 

I.m 

The entropy is in turn given by 

S = -A^E'TrlMn/k + ( l ~ /OMl - /0] , (38) 

k 

where /k is the matrix of expectation values defined in Eq. (|35| ) . 
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B. Zero temperature 

At zero temperature, the variational state gives an approximate ground state \G) and 
the thermal expectation values evolve into expectation values with respect to this ground 
state i.e., ((\l/ k <8> ^ k )) * s replaced by (G]^ ® \l/ k |Cr). The zero-temperature self-consistent 
equations are obtained from Eq. (^l|) in the limit /3 — > oo. 

Combining the facts that the ground state is the vacuum state for the quasiparticles and 
that the first (last) four elements of x k are annihilation (creation) operators, we can identify 
the diagonal matrix g = (G|(xl ® Xk)\G)- Its diagonal entries are (0, 0, 0, 0, 1, 1, 1, 1), and 
from Eq. ( |3"4"D with (3 — > oo it follows that four of the eigenvalues of h(k) must be negative 
and the rest positive. We then solve for the unitary matrices ?7 k that diagonalize the /i(k)'s 
defined in Eq. (E9I), in such a way that the eigenvalues in the diagonal matrices _D k are in 
descending order. The self-consistent equations to be iterated are then Eq. (p9| ) and 

A - = lj2"^MBmUt9U k ) , (39a) 

8 k 

D k = U k h(k)Ul (39b) 



C. An alternative set of equations for the ground state at half filling 

The self-consistent equations must be solved numerically. Although the equations are 
formally simple, it is quite challenging to numerically carry out a search of the solution 
space. We therefore present an alternative method to find the ground state, and use it to 
derive a theorem of stability of the superconducting and the antiferromagnetic solutions at 
half filling for some regimes of U and J. 

Naively, the expectation value of the Hamiltonian, Eq. (|19|), is a simple quadratic form 
and one should just find its minimum. However, since the expectation values (a^ 1 ) are 
constrained by the fact that they represent expectation values of fermion operators, the 
terms cannot be independently varied, and there are constraints on the set of (a^). These 
restrictions were automatically satisfied in the previous analysis, since a k was explicitly 
computed through a canonical transformation. Another way to proceed is to try to find a 
constrained quadratic minimum directly without first calculating a canonical transformation. 

To express the fermionic constraints, we define the following matrix 

A k = £ (a£)(aZ'){B m ,B m ,}, (40) 

where {A, B} denotes the anticommutator {A, B} = AB + BA. The matrix Bq = B® Q = — II 
that is excluded from the sum is the only basis matrix with nonzero trace. First we state 
the lemma that reformulates the problem of minimizing the energy. 

Lemma 1 The restrictions on the expectation values (a.^) for the variational solutions of 
Eq. are: E m «) 2 = \, (Ko)k> = ~\, and A k = V k. 
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This lemma is proved in App. ||. The importance of the lemma is that it shows that the 
minimization problem of Eq. (19) is a quadratic minimum subject to quadratic constraints. 



This enables us to search for minima of the unconstrained problem, which, if they are found 
to satisfy the constraints, must also be minima of the constrained problem. A class of such 
solutions are introduced in the following theorem (proved in App. |B|) and corollary. 

Theorem 1 Consider the Hamiltonian in Eq. flTQ ) with fixed coefficients a^ m and b^ m . De- 
fine M as the subset of purely quadratic and irreducible terms, M = {m : a^ m = VZ , b^ m = 
b m b~i,i(m)} ! where l(m) is a function which attaches one single I to each m. Assume there ex- 
ists anmeU, such that (b^ < and ^l/^W^I < &m|/£ M /# m) | Vk, k' ; Vm G M\{m}). 
If ({a™) = 0Vk, Vm G A/"\ fa) is a sufficient condition for (A^ = OVkX then the same 
(a) 's will be zero also in the minimizing solution of the constrained problem. 

In the case of half-filling the following corollary now follows for the two important cases of 
antiferromagnetic and s-wave superconducting ordering. 

Corollary 1 The state of lowest energy at half filling (\x = 0) for the Hamiltonian (H = 
H + H Uuhh + H Ucis ) in Eqs. flT^), (ff/p and (|7^ is AF ((a%) ^ o) ifO < J <U, and s-wave 
SC-CDW ((a 2 J ■£ 0) if U < J < -U/3. 

The theorem follows by making two observations. First, the possible low-energy states follow 
from Theorem [I] by inspection of the prefactors of the quadratic terms in the Hamiltonian, 
and by verifying that A^ = if all (a)'s are zero except the hopping (a^o) anc ^ either 
(af ) or (otot)- Secondly, we observe that the problem that results by setting all other order 
parameters to zero is analogous to the ordinary s-wave SC case, where it is well known that 
any attractive interaction results in a finite order parameter. 

We further note that there are two degenerate superconductivity solutions if the con- 
straints are disregarded. One is the ordinary s-wave superconductor, and the other is the 
staggered superconductor. However, the latter is ruled out by the fact that the constraint 
Ak = is not fulfilled. Apart from predicting these low energy states, Theorem [l] also proves 
the stability of these phases with respect to small perturbations to the Hamiltonian. Since 
the ferromagnetic state has A^ ^ 0, ferromagnetic ordering and nonzero hopping cannot be 
present simultaneously, at least not at the same location in k-space. The predictions of the 
Corollary are illustrated in Fig. [| 

VII. NUMERICAL METHODS 

To find the phase at a (U, J, fi, T)-point in a phase diagram, we choose an initial value of 
the pitch vector q, solve the self-consistent equations Eqs. (|29|), (|31|) and (|35|) numerically, 
and calculate F(q) using Eqs. (|37D and (|38|) . 

The pitch vector q is then varied to find which q gives the solution to the self-consistent 
equations with the lowest value of F(q). The non-zero mean-field order parameters A; im 
defines, together with q, the particular phase for the (U, J, /i, T)-point in parameter space. 

Complete phase diagrams are obtained by the following two steps: 

• The parameters (U, J, /z, T, q) are swept, with an initial set of A m 's generated at ran- 
dom and the self-consistent equations are iterated until a fixed point is reached. This 
gives us a rough picture of the states that are present in the phase diagram. 
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• The accuracy of the boundaries between the phases in the phase diagram is improved. 
Here the self-consistent equations are solved using Broyden's, methodic which is often 
more efficient than the previous iterative method. 

To find Ai >m for a particular value of q, we cover the reduced Brillouin zone by a discrete 
lattice. Care has to be taken not to break any of the symmetries of the problem. Fig. |T| 
shows how a 32-point lattice is laid out. The most time-consuming numerical step is to 
diagonalize the 8x8 matrix in the argument of the exponential function in Eq. ( pop at every 
k-point; this must be done for each iteration. 

Choosing random order parameters as initial conditions for recursion is useful when there 
is no a priori information about the expected solutions of the self-consistent equations. A 
complication of this method, however, is that the iteration tends to fall into cycles. We cured 
this by including a tail of exponentially damped previous iterates at each step. Sometimes, 
the procedure still did not converge, and several initial points must be used before a fixed 
point was found. 

To obtain a complete \i-T phase diagram for fixed U and J, we cover the (/i, T, q) space 
with roughly 1500 points, and repeat the iterative procedure 10 times. On an ordinary 
workstation it takes of the order of a week of CPU time to trace out the phase diagram 
using 98 points in the reduced Brillouin zone and solving the self-consistent equations to an 
accuracy of 1 %. Of course, solutions could be missed by chance since we use random initial 
guesses, and phases occurring in narrow regions of the phase space could be missed since 
the parameter space is not covered with a fine enough mesh. 

After the different phases have been identified, we obtain more accurate solutions of 
the self-consistent equations using Broyden's method. This method cannot be used from 
the beginning since the initial guess has to be close to the final answer for the method to 
converge. The method is also slow if the number of order parameters is very large. Here, we 
therefore eliminate from the Hamiltonian all order parameters that are known to be zero in 
the corresponding regions. 

If the state of lowest free energy has q = the minimizing solution can be obtained 
directly by Broyden's method, and in this case the phase boundaries are located to high 
accuracy; generally by using 800 points in the reduced Brillouin zone. 

If q is non-zero, there is the extra problem of minimizing the free energy with respect to q. 
We do this by using Broyden's method for solving the self-consistent equations and extending 
it with a simple numerical algorithm that also minimizes in q using Brent's methodic To 
obtain results within reasonable computing time, most of the spiral spin-wave calculations 
have been performed using a 392-point lattice in the reduced Brillouin zone. 

In order to allow for the most general spiral solutions, no restrictions should be imposed 
on the spiral spin-wave parameter q, but that would make the problem numerically un- 
manageable. Instead we have focused on the question of whether the low-energy state is a 
spiral spin- wave or not. Assuming the spiral spin- wave not to break the lattice symmetries 
completely, the quantization axis q could be twisted either in the (1, 1) or (1,0) direction. 
We further concentrated on the latter case case, q = (q, 0), since it turned out to give a 
slightly lower free energy than the diagonal twist in some regions of the phase diagram. 
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VIII. PHASE DIAGRAMS 

In order to present a set of complete phase diagrams for the extended Hubbard model, 
we would have to probe all combinations of values of the four parameters U, J, filling (or 
Id) and T. This is an unfeasible task, and we restrict ourselves to certain cross sections that 
we hope to capture generic behavior. Due to the particle-hole symmetry as /i — > —fi, we 
restrict our phase diagrams to hole doping (/i < 0). We can further set t = 1 without loss 
of generality. 

A. Zero temperature and half filling 

The most fundamental cross section is (T = 0, ji = 0) corresponding to the ground state 
of the half-filled extended Hubbard model in Fig. |3|. In this case, when U and J are both 
positive, as well as for a region of negative U, the antiferromagnetic (AF) state forms a 
numerically very stable solution. This is consistent with Corollary |1|. The corollary further 
indicates a region (SC-CDW) of degenerate s-wave SC and charge- density- wave state, which 
is confirmed by the numerical simulations for negative U and intermediate J. 

An interesting feature occurring in a region outside the validity of the assumptions of 
Corollary [TJ, but which is numerically very robust, is the tongue of (i-wave SC-CDW, (s + id)- 
wave SC-CDW and ((i-wave SC-CDW) +AF between the AF and the s-wave states. The 
d-w&ve order parameter is dominant along the center line of the tongue while it vanishes at 
the boundaries. Numerically we also see that the s and d-w&ve order parameters form parallel 
pseudospin vectors, and the explicit form of the base matrices of these vectors indicate that 
the mixed state breaks time reversal symmetry, i.e., it is an s + id state. This is consistent 
with the constraint A^ = in Lemma [l], since the constraint requires the two pseudospin 
vectors to be parallel unless some other order parameter is non-zero. 

The ferromagnetic (FM) zone that is seen for negative J is just barely numerically stable 
close to the phase boundaries. It is, however, quite stable deeper inside the FM region. The 
numerical difficulties near the FM phase boundaries can be understood by the constraint 
conditions discussed in Lemma [TJ. Since, in the absence of a third order parameter, hopping 
and FM order cannot exist simultaneously at the same point in the Brillouin zone, there will 
be distinct regions in k-space. The FM ordering occurs close to the Fermi surface, while the 
hopping expectation value is finite near the origin in reciprocal space. The sharp boundary 
between the two regions results in discontinuities in the numerical solution. 

B. Finite temperature 

We now turn to phase diagrams at non-zero temperature and variable filling for some 
particular fixed values of U and J. Estimates using physical models for the high-T c materials 
have suggested that U should be of the order 5 with J much smaller. The relevance of such 
estimates for the present calculation is questionable since they were made for models of 
high-T c materials that include other types of interactions such as nearest-neighbor charge 
repulsion. Moreover, a large fraction of our Heisenberg term could be considered as coming 
from effective second-order corrections to the Hubbard interaction, which are otherwise 
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neglected in our mean-field approach. Therefore, we concentrate on parameter values that 
give interesting phenomena. 

There are two features that we are particularly interested in exploring. The first is d- 
wave superconductivity, and the second is spiral spin-waves. To keep some connection to 
high-T c materials, we require the model to be antiferromagnetic at half filling, and therefore 
U and J should be positive. In order to study <i-wave superconductivity, U should not be 
too large. The spiral spin-wave states, on the other hand, have been observed for the pure 
Hubbard model with intermediate U, and in this case J should be zero or at least small. 

We start out by investigating the phase diagram for U = 0.5 and J = 2 which contains 
a zone of <i-wave superconductivity (see Fig. [|). The system is antiferromagnetic close 
to half filling. This state persists up to the Neel temperature, where there is a second- 
order phase transition to the normal state (NS) that has no broken symmetry. At low 
temperatures and moderate doping, we have a d-wave superconductor which is separated 
from the antiferromagnet by a first-order phase transition. For temperatures higher than the 
SC critical temperature, there is a first-order phase transition between the antiferromagnet 
and the normal state. This first-order boundary terminates at higher temperatures at a 
critical point, where the transition becomes second order. No qualitative differences between 
this phase diagram and that for U — is observed. 

Both antiferromagnetic phases and s- and d-wave superconductors have been observed in 
other calculations.E-a. Those studies did not consider the possibility that the phase transition 
could be first order, but by carefully comparing the free energies we have observed this type 
of transition between the AF and the <i-wave SC states at low temperature. If, as in our 
model, the total electron number is fixed, the first-order phase transition results in phase 
separation. 

By increasing J and U, the antiferromagnetic region is enhanced. The ci-wave zone 
grows with increasing J, but is also shifted to larger doping as the antiferromagnetic region 
expands at the same time. The overall size of the SC region in the T — jjl phase diagram is 
insensitive to a change in U. The s-wave superconducting zone is also enhanced by larger 
J, while it diminishes if U is increased too much. It should also be noted that the AF and 
the s-wave SC regions gain more from an increase in J than the d-w&ve regions do. This 
is illustrated in Fig. || where the phase diagram for U = and J = 4 is exhibited. If J is 
sufficiently large, we expect the <i-wave to disappear in favor of the s-wave SC and the AF 
phases. 

For large U and small J, spiral spin- waves appear as shown in Fig. [] for J = 0. The 
phase diagram is shown in Fig. ||, where SSW indicates the spiral spin-wave with the pitch 
(n — q, 7r), which is obtained by applying the twist q = (1, 0)q to the AF (ir, ir) state. At low 
temperature, the SSW state is separated from the AF state by a first-order phase transition 
with a wide coexistence region as a function of density. For more elevated temperatures, 
the separation line becomes second-order and here the spiral pitch parameter q goes to zero 
as the phase boundary is approached. On the contrary, along the phase transition from the 
SSW to the normal state the spiral magnetic order parameter vanishes in magnitude while 
q stays finite. No superconductivity is seen in Fig. ||. The s-wave superconducting state is 
suppressed by the large value of U, and the <i-wave state is not seen either since J is zero. 

The energy gain of the spiral spin- waves is numerically very small, of the order of a 
hundredth of the overall condensation energy. The introduction of a Heisenberg interaction 
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could therefore have a large influence. To investigate this issue, we introduce a small J = 0.1 
while keeping U = 5, and compare Fig. |6| with Fig. |7|. What was a second order transition 
between the SSW and the normal state for J = has now become first order. There is also 
a temperature range at which the AF is "reentrant" as a function of doping, and where the 
AF-SSW phase boundary move toward lower temperatures and becomes first order. When 
increasing J, this line of phase transitions rapidly migrate towards lower doping and very 
soon the whole spiral spin-wave region is gone. This explains why no spiral spin-wave is 
seen in the phase diagrams for small U and large J. 

For positive J and sufficiently low densities, the antiferromagnet and the spiral spin-wave 
must eventually disappear, leaving room for the s-wave and d-wave superconducting states 
which may persist the rest of the way to zero filling. The critical temperature decreases 
rapidly with decreasing J and we were unable to confirm this scenario numerically. 

IX. DISCUSSION 

The Hubbard model serves as a simple model for high-T c superconductors. Although its 
simple appearance, the model is still poorly understood and many sophisticated techniques 
for studying specific features of the model have been proposed in the literature. As a guide 
to this realm of possibilities, it is important to have a good understanding of all possibilities 
that a "simple" mean-field analysis can provide. We have therefore used a generalized 
Hartree-Fock-Bogoliubov theory and numerical simulations to compute phase diagrams for 
the extended Hubbard model. All the conventional order parameters, like s and <i-wave 
superconductivity, charge-density waves, and Neel and spiral antiferromagnetic states states, 
have been included in one unifying framework, making no a priori assumptions about the 
nature of the broken symmetries. We have further shown that, in mean-field theory, no new 
mixed phases arise at finite doping and temperature in the extended Hubbard model with 
positive values of U and J. In our investigation, we have seen the time-reversal symmetry- 
breaking superconducting phase s + id only in a narrow region with negative U. Close to 
this region there is also a region of mixed antiferromagnetism and d-w&ve superconductivity. 

Our method allow phase separation to occur, which it also does in certain regions. The 
energy differences that we find between Neel and spiral antiferromagnets is so small, that we 
do not want to make any strong statements about whether phase coexistence would survive 
a more refined analysis or not. However, the energy difference between the AF and the 
normal state at the first order phase boundary is substantial. A phase separation between 
these two states has also been suggested both from theoretical and experimental groundsill 
Another thing to keep in mind is that we require the total number of electrons to be fixed, 
while in the high-T c materials there are large charge reservoirs surrounding the 2D planes 
that are perhaps better modeled by a fixed chemical potential. Under such circumstances, 
the phase coexistence may well be suppressed. 

There have been several earlier studies of spiral spin-waves for the Hubbard model ex- 
ploiting slave boson and ordinary Hartree-Fock techniques. Most of these studies have 
concentrated on zero temperature, and our corresponding results are consistent with those. 
However, we have also extended the analysis to finite temperature. 

Another new phenomenon that we have studied is how the spiral spin-waves are affected 
by a nearest-neighbor Heisenberg term in the Hamiltonian. We have observed that the (1, 0) 
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spiral spin- wave phase is easily destroyed upon the introduction of a positive- J Heisenberg 
interaction. We have mostly concentrated on spin-waves in the (1, 0)-direction since we 
found that the energy differences are very insensitive to the pitch direction, at least for 
the pure Hubbard model. To make the study complete, other spin directions should also 
be studied. However, it is likely that these small energy differences are insignificant with 
respect to the overall crudeness of our analysis, although one might argue that their relative 
difference is to be taken seriously. 
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APPENDIX A: THE POINT GROUP C AV 

The point-group symmetry of the 2D square lattice is C Av , since the lattice is invariant 
under 90° rotations around the z-axes and under reflections in the lines v and v ' in Fig. || 
The group elements of C iv are the identity (/), 90° rotations (C 4 ), 180° rotations (Cf), 
reflections in v (a v ) and reflections in v' (<v). This group has four ID irreps (A\, A 2 , £>i, 
B2) and one 2D irrep E. The character table together with examples of basis functions for 
the different irreps are given in Table O. Our main use of the C± v irreps is to distinguish 
between d and s-wave superconductivity order parameters. The s-wave ordering has the full 
symmetry of the lattice, i. e., it belongs to the Ai representation. The <i-wave ordering, on the 
other hand, is antisymmetric under reflections in v ', and belongs to the B\ representation. If 
we would see any p-wave states, these would belong to the 2D E representation since these 
states are antisymmetric under parity ((x,y) —* (—x, —y))- 



APPENDIX B: PROOFS OF THE ZERO-TEMPERATURE, HALF-FILLING 

THEOREMS 

In this appendix we give the proofs of the theorems concerning the minimization of the 
ground-state energy. A bunch of related theorems have also been derived by Bach et «/.e3 
First we prove the lemma for how the energy minimization problem can be recast into the 
problem of minimizing the expectation value of the Hamiltonian written in terms of (at^) . 

Lemma |l] The restrictions on the expectation values (ac™) for the variational solutions of 
Eq. QI%) are: E m «) 2 = \, («>)k> = ~\, ^d A k = V k. 

Proof. The space of all possible variational solutions is defined by the constraint that the 
Hermitian matrix (G|^k <8> ^kl^) nas t ne four-fold degenerate eigenvalues and 1, since 
it has the same eigenvalues as g = (G\(xL ® Xk)l^)- O ur a i m i s to find the corresponding 
constraints on the coefficients (a k ) in the expansion 

,t 



(G|*L®*k|G> = £«>2Jm- (Bl) 



First of all, since B$ = Bq = — 11 is the only basis matrix with a non-zero trace, one has 
((. a oo)k) = —\- Let us next define the trace-less matrix X k 

X k = (G|*t®tt k |G)-§l. (B2) 

This matrix has the fourfold degenerate eigenvalues ±| and the expansion 

X k =£«)£ m . (B3) 

Furthermore, {X k ,X k } has the 8-fold degenerate eigenvalue |, meaning that {X k ,X k } = 
^11, so that using Eq. (|9]) we have 

±l = {X k ,X k } = ^2«) 2 lL + ,4 k , (B4) 

where v4 k is the non-diagonal part defined in Eq. (flOf) . Since A k is orthogonal to 11, it follows 
that v4 k = and that 

X^/^™\ 2 = /n,9\ 2 4- >T /™. m \ 2 = r-I^ 2 4-1 = 1 



Z^\ a k 



) 2 = K) 2 +E«) 2 = H) 2 + i = l- (B5) 



m^O 



Q. E. D. 

For some specific solutions, this lemma leads to the following theorem which considerably 
simplifies the search for solutions. Here we consider the larger space of solutions that arises 
if we disregard all constraints except the normalization J2 m ( a k) 2 — §• ^ solution of the new 
problem that happens to fulfill all the constraints, must then be a solution of the original 
problem as well. 

Theorem |l| Consider the Hamiltonian in Eq. (Tl) with fixed coefficients a^ m and b^ m . De- 
fine M as the subset of purely quadratic and irreducible terms, N = {m : a^ m = VZ , b^ m = 
b m 8i,i(m)}, where l(m) is a function which attaches one single I to each m. Assume there ex- 
ists anm G M, such that fa < andb^ 1 ^ ^ ] \ < 6 m |/3 k M /3 k ( , m) | Vk,k', Vm G M\{m}). 
If ({a™) = 0Vk, Vm E Af \ rh) is a sufficient condition for (A k = OVkj, then the same 
(a) 's will be zero also in the minimizing solution of the constrained problem. 

Proof. The lowest-energy solution should be minimal under any variation of (a)'s that fulfills 
the normalization constraint (|B5|). Since there is only one constraint per k in the simplified 
problem, it is possible to keep all (a)'s fixed except two and still fulfill the constraint. Sup- 
pose that we vary only (a^) and (« k ), where m,n G A/". A variation around a minimizing 
solution must then fulfill 

+bj^& n \al,)5(al) = 0, (B6) 

(<><K<> + <<>%£> = 0. 

From our assumptions, we have bfh < 0, and to consider competing solutions one must also 
have b n < 0. Since there are no constraints imposed on the signs of (a™) and (a k ), it 
is obvious that a minimizing solution is obtained by choosing sgn((a k a )) = sgn(/5 k ) and 
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similarly for (af k ). Taking these sign considerations into account and eliminating S{a^) and 
5(a£) yields 

E"(^l/5r ) ^ ) ll(^)K)l 

-&„|rf n) ^ n) ||K,)(«ni) = 0. (B7) 

Summing this equation over k gives 

£>^fW } | - bnl^^Dla^l = , (B8) 

k,k' 

and since from our assumption, &m|/3 k /3 k / | < & n |/3 k /3 k / |Vk, k', the solution must be 
either a^ = OVk or «£ = OVk. Of these two, a k = is obviously the solution of lowest 
energy. This argument is then applied to all purely quadratic terms n G J\f\m. 
Q. E. D. 
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FIGURES 
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-71 k x * 

FIG. 1. The first Brillouin zone of the square lattice with lattice constant 1. The line ek = 
is indicated as well as the k-points in the reduced zone that are used in a numerical simulation 



with 32 k-points (Sec. VII). The contributions from the points on the /c x -axes are scaled down by 
a factor of two since they would otherwise be overcounted. 




FIG. 2. The low energy states of the extended Hubbard model according to Corollary [l] is 
shown in grey. The white areas are indeterminate from the corollary, and the phases here must be 
computed numerically. 
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FIG. 3. Phase diagram at half filling and zero temperature as a function of J and U. Sec- 
ond-order phase boundaries are drawn as full lines, while first order phase boundaries are drawn 
as dashed lines. The dotted lines are extrapolations of the numerically derived full (dashed)-line 
phase boundaries. 
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FIG. 4. Phase diagram for U = 0.5, J = 2. Second-order phase transitions are drawn as full 
lines, while the boundaries of regions with phase separation are drawn as dashed lines. The regions 
of phase separation are denoted A / B where A and B are the two coexisting phases. Regions of 
s-wave and d-wave SC are indicated by s-w and d-w, respectively. 
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FIG. 5. Phase diagram for U = and J = 4. The d-wave and s-wave superconducting regions 
are separated by a first order phase transition, with a very narrow coexistence region. The inset 
shows a magnification of the coexistence region. 
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FIG. 6. Phase diagram for U = 5 and J = 0. The spiral spin-wave with pitch (tt — q, n) is 
denoted SSW. 
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FIG. 8. The symmetry axes v and v' of the square lattice. 



26 



TABLES 



TABLE I. Construction of the 8 x 8 (3a matrices from the 4x4 Dirac matrices expressed in 
terms of Pauli matrices. The index j runs from 1 to 3. The notation 7T denotes complex conjugate 
(not adjoint). 



Pauli 

4x4 
Dirac 

8x8 

Pa 



{0-1,0-2,0-3) 



To 



11 
-1 



7o 






-To 



(3 j+3 = iZfofyZ 



1 

1 



7j = 



-i 

i 

Oj 
-oj 

n °* 

7* 



1 
-1 



TABLE II. The character table of the point group C± v together with examples of basis functions 
for the different irreducible representations. 





I 


cl 


C A 


a v 


a v i 


Examples of functions 


Ai 


1 


1 


1 


1 


1 


1, cosk x + cos k y 


A 2 


1 


1 


1 


-1 


-1 


sin 2k x sin k y — sin 2k y sin k x 


Bi 


1 


1 


-1 


1 


-1 


cos k x — cos k y 


B 2 


1 


1 


-1 


-1 


1 


sin k x sin k y 


E 


2 


-2 











{smk x , smky} 



27 



